## Code that creates figure presenting results by crisis severity. 
## Last changed: 2025-05-23

sink("C:/Userdata/Shared/Dofiles/DoAnalysis/JanKalle/ReplicationFiles/AnalysisScripts/Output/FigureA19.txt", split=TRUE)


library(PSweight)
library(cobalt)
library(MatchIt)
library(distances)
library(moreparty)
library(party)
library(caret)
library(bonsai)
library(randomForest)
library(tidymodels)
library(tidyverse)
library(rmarkdown)
library(tinytex)
library(tibble)
library(knitr)
library(car)
library(fst)
library(stargazer)
library(pryr)
library(stringr)
library(fst)
library(lubridate)
library(data.table)
library(summarytools)
library(ggplot2)
library(cowplot)
library(ggpubr)
library(haven)
library(broom)
library(fixest)
library(modelsummary)
library(gt)
library(MASS)
library(survey)
library(collapse)
library(readxl)
library(Hmisc)
library(DescTools)
library(ranger)
library(stablelearner)
library(future.apply)
library(lme4)
library(partykit)
library(ggiplot)

  LA91 <- as.data.table(read_xls("E:/ProjData/Jan_Kalle/la1991.xls"))
  LA91[, Kommun := as.numeric(KommunKod)]
  
  RTB94 <- read.fst("D:/SCB_ConPol/Fst/RTB/RTB_1994.fst", as.data.table=T)
  RTB94 <- RTB94[, .(LopNr, Kommun)]

  
  RTB91 <- read.fst("D:/SCB_ConPol/Fst/RTB/RTB_1991.fst", as.data.table=T)
  RTB91 <- RTB91[, .(LopNr, Kommun)]
  
  RTB91 <- LA91[, .(LAKod, Kommun)][RTB91, on = "Kommun"]
  RTB94 <- LA91[, .(LAKod, Kommun)][RTB94, on = "Kommun"]
  
  RTB94[Kommun == 1535, LAKod := 46] # Bollebygd - Borås LA
  RTB94[Kommun == 1814, LAKod := 64] # Lekeberg - Örebro LA
  
  FodUppg <- read.fst("D:/SCB_ConPol/Fst/RTB/FodelseUppg_2024.fst", as.data.table=T)
  FodUppg <- FodUppg[, cID := .N, by = LopNr][cID == 1 & !is.na(LopNr), ][, cID := NULL]
  FodUppg <- FodUppg[, .(LopNr, FodArMan, UtlSvBakg, Kon)]
  FodUppg[, FodAr := floor(as.numeric(FodArMan)/100)]
  FodUppg[, FodArMan := NULL]

  #deso91 <- read_fst(paste0("D:/SCB_ConPol/Fst/SAMS/DESO_", 1991, ".fst"), as.data.table=T)
  #deso91 <- deso91[, cID := .N, by = LopNr][cID == 1 & !is.na(LopNr), ][, cID := NULL]
  #deso91[, Kommun := str_sub(Deso, 1, 4)]
  
  #deso94 <- read_fst(paste0("D:/SCB_ConPol/Fst/SAMS/DESO_", 1994, ".fst"), as.data.table=T)
  #deso94 <- deso94[, cID := .N, by = LopNr][cID == 1 & !is.na(LopNr), ][, cID := NULL]
  #deso94[, Kommun := str_sub(Deso, 1, 4)]
  
  L91 <- read_fst(paste0("D:/SCB_ConPol/Fst/Lisa/Lisa_1991.fst"), as.data.table = T)
  L91[, cID := .N, by = lopnr][cID == 1 & !is.na(lopnr), ][, cID := NULL]
  setnames(L91, "lopnr", "LopNr")
  L91[, Ar := 1991]
  L91[, Sysselsatt := as.numeric(SyssStatG == 1)]
  #L91 <- deso91[, .(LopNr, Kommun)][L91, on = "LopNr"]
  L91 <- RTB91[, .(LopNr, LAKod)][L91, on = "LopNr"]
  L91 <- FodUppg[L91, on = "LopNr"]
  
  L94 <- read_fst(paste0("D:/SCB_ConPol/Fst/Lisa/Lisa_1994.fst"), as.data.table = T)
  L94[, cID := .N, by = lopnr][cID == 1 & !is.na(lopnr), ][, cID := NULL]
  setnames(L94, "lopnr", "LopNr")
  setnames(L94, "SyssStat", "SyssStatG")
  L94[, Ar := 1994]
  L94[, Sysselsatt := as.numeric(SyssStatG == 1)]
  L94[is.na(Sysselsatt), Sysselsatt := 0]
  L94 <- RTB94[, .(LopNr, LAKod)][L94, on = "LopNr"]
  #L94 <- deso94[, .(LopNr, Kommun)][L94, on = "LopNr"]
  L94 <- FodUppg[L94, on = "LopNr"]
  
  L91 <- L91[, .(LopNr, LAKod, Sysselsatt, FodAr)]
  L94 <- L94[, .(LopNr, LAKod, Sysselsatt, FodAr)]
  L91[, Alder := 1991-FodAr]
  L94[, Alder := 1994-FodAr]
  
  L91[between(Alder, 25, 64) & !is.na(LAKod), LASyss91 := mean(Sysselsatt, na.rm=T), by=LAKod ]
  L94[between(Alder, 25, 64) & !is.na(LAKod), LASyss94 := mean(Sysselsatt, na.rm=T), by=LAKod ]
  Kom91 <- unique(L91[!is.na(LAKod) & !is.na(LASyss91), .(LAKod, LASyss91)])
  Kom94 <- unique(L94[!is.na(LAKod) & !is.na(LASyss94), .(LAKod, LASyss94)])
  
  LASyss9194 <- Kom91[Kom94, on = "LAKod"]
  LASyss9194[, DSyss9491 := LASyss94-LASyss91]
  qsu(LASyss9194)
  
  RTB91 <- LASyss9194[RTB91, on = "LAKod"]
  RTB91 <- FodUppg[RTB91, on = "LopNr"]
  kLA <- quantile(RTB91[between(FodAr, 1942, 1967), DSyss9491], na.rm=T)
  
  RTB91[between(DSyss9491, kLA[1], kLA[2]), KrisLA := 2]
  RTB91[between(DSyss9491, kLA[2], kLA[3]), KrisLA := 2]
  RTB91[between(DSyss9491, kLA[3], kLA[4]), KrisLA := 1]
  RTB91[between(DSyss9491, kLA[4], kLA[5]), KrisLA := 1]
  
  # Specify the name of the treatment variable
  tvar <- "T94_90"
  
  # Specify the the birth cohorts to be used for the analysis
  lbyr <- 1942
  ubyr <- 1967
  
  # Specify if only include Swedish citizens 
  citizen <- 1
  
  # Specify the propensity score specification
  psformula <- formula(T94 ~
                         Kon + DAStod_1982 + DAStod_1983 + DAStod_1984 +
                         DAStod_1985 + DAStod_1986 + DAStod_1987 + DAStod_1988 +
                         DAStod_1989 +
                         DAntBarn_91 + Sambo_91 + InvBak  |
                         SSYK3_90 + SNI2_91 + Kommun_91 +
                         rCARB_1982 + rCARB_1983 + rCARB_1984 + rCARB_1985 +
                         rCARB_1986 + rCARB_1987 + rCARB_1988 + rCARB_1989 +
                         pArbInk_1990 + pArbInk_1991 + UtbAr_91 +
                         FodAr + Nom82 + Nom85 + Nom88 + Nom91)
  
  ## Estimate two models, one for the complete sample including all observations,
  ## the second on the restricted sample only including individuals for which
  ## we can observe all elections bw 82-18 in which they are eligible to run for office. 
  
  Models <-  vector('list', 10)
  j <- 1
  urval <- 1
  
  for(kk in c(1, 2)) { 
    for(i in c(0, 4)) { 
      # Read in the data sets
      NomVald <- read_fst("E:/ProjData/Jan_Kalle/JoPReplications/DB_NomVald.fst", as.data.table=T)
      RTB <- read_fst("E:/ProjData/Jan_Kalle/JoPReplications/DB_RTB.fst", as.data.table=T)
      db <- read_fst("E:/ProjData/Jan_Kalle/JoPReplications/DB_Nom.fst", as.data.table=T)
      
      db[, T94 := get(tvar)]
      mdb <- db[between(FodAr, lbyr, ubyr) & (SvMedb82_18 >= citizen) & (AllaVal >= urval),]
      rm(db)
      gc()
      
      
      # Split the in SES quartiles
      mdb[!is.na(avSES3), ravSES := cut_number(avSES3, 4, labels = FALSE)]
      if(i==1) mdb <- mdb[between(ravSES, 1, 1)]
      if(i==2) mdb <- mdb[between(ravSES, 2, 2)]
      if(i==3) mdb <- mdb[between(ravSES, 3, 3)]
      if(i==4) mdb <- mdb[between(ravSES, 4, 4)]
      
      mdb <- RTB91[, .(LopNr, KrisLA)][mdb, on = "LopNr"]
      mdb <- mdb[KrisLA == kk]
      
      # Use logit model to estimate propensity scores
      mdb <- mdb[complete.cases(mdb),]
      tm <- proc.time()
      logmod <- feglm(psformula,  
                      data = mdb, family="binomial")
      proc.time()-tm
      summary(logmod)
      
      # Extract predicted probability (propensity scores)
      if(logmod$nobs == logmod$nobs_origin){ 
        mdb[, pT94 := predict(logmod, type = "response")]
      }else{
        mdb[logmod$obs_selection$obsRemoved, pT94 := predict(logmod, type = "response")]
      }
      
      mdb <- mdb[complete.cases(mdb),]
      
      # Now use the Full matching approach of Sävje et al. to obtain ATT and ATE weights.  
      
      tmp <- proc.time()
      set.seed(7892)
      m.out <- matchit(T94~1, data = mdb,
                       distance = mdb$pT94, method = "quick", estimand = "ATT")
      proc.time()-tmp
      m.data1 <- as.data.table(match.data(m.out))[, .(LopNr, weights)]
      setnames(m.data1, "weights", "gm_att")
      
      mdb <- m.data1[, .(LopNr, gm_att)][mdb, on = "LopNr"]
      gc()
      RTB <- mdb[RTB, on = "LopNr"]
      
      rm(mdb)
      gc()
      
      rm(list=setdiff(ls(), c("NomVald", "RTB", "rtbvar", "Models", "j", "tvar", "lbyr", "ubyr", "citizen", "psformula", "urval", "RTB91", "kk", "i")))
      gc()
      
      
      mdb <- RTB[!is.na(gm_att), 
                 .(Kon, Nom82, Nom85, Nom88, Nom91, Sysselsatt_85, DAntBarn_91, Sambo_91, 
                   DStud_91, DForLed_91, DSjukRe_91, DSocBidrFam_91, DBostBidrFam_91,
                   UtbAr_91, FodAr, InvBak, isei, Yrke80_90, SNI2_90, Kommun_91, fp_LoneInk_85, fp_LoneInk_90, fp_LoneInk_91,
                   gm_att, Nom, LopNr, T94, Ar, pT94)]
      
      gc()
      
      mdb[, Nom := Nom*100]
      mdb[, cID := .N, by=.(LopNr, Ar)]
      #mdb <- mdb[(Ar-FodAr)<70, ]
      
      rm(RTB)
      gc()
      
      
      Models[[j]] <- feols(Nom~T94+i(Ar, T94, ref=1991) | Ar, cluster ="LopNr", 
                           weights = mdb[, gm_att], 
                           mdb[])
      
      
      print(qsu(mdb[Ar==1991, .(Nom)]))
      print(summary(Models[[j]]))
      print(uniqueN(mdb[T94==0, LopNr]))
      print(uniqueN(mdb[T94==1, LopNr]))
      print(uniqueN(mdb[, LopNr]))
      print(dim(mdb))
      print(qsu(mdb$cID))
      
      print(i)
      print(kk)
      j <- j+1
      
    }  
  }  
   
  valar <- c(1982, 1985, 1988, 1991, 1994, 1998, 2002, 2006, 2010, 2014, 2018, 2022)
  options(scipen=999)
  
    for(i in 1:4) {
      
      ses <- ""
      if(i %in% c(1)) ses <- paste0("High Crisis, All SES") 
      if(i %in% c(2)) ses <- paste0("High Crisis, SES Q4")
      if(i %in% c(3)) ses <- paste0("Low Crisis, All SES") 
      if(i %in% c(4)) ses <- paste0("Low Crisis, SES Q4")
      
      gi <- ggiplot(Models[[i]]) 
      gp <- ggplot(gi$data, aes(x, estimate)) +
        geom_point() + 
        geom_hline(yintercept=0, linetype="longdash") +
        geom_errorbar(aes(ymin=ci_low,ymax=ci_high), width=0.4, size=.6) +
        geom_line(linetype="dotted") +
        scale_x_continuous(breaks = valar,  labels = valar, limits = c(1981, 2023)) +
        scale_y_continuous(limits = c(-.5, 1.15), breaks = round(seq(-.5, 1.1, by = .2), digits=1)) +
        xlab("Election Year") +
        ylab("Effect on Candidacy") +
        ggtitle(ses) +
        theme_classic(base_size = 14) +
        theme(legend.title=element_blank()) +
        theme(legend.position = "bottom", panel.grid.major.y = element_line(color = "grey95"))
      
      gp
      #ggsave(paste0("C:/Userdata/Shared/Dofiles/DoAnalysis/JanKalle/Figures/Fig_", j, "_", i, "kk.pdf"), width = 10, height =7.5, units = "in") 
      
      assign(paste0("gp", i), gp)
    }  
  
  ggarrange(gp1, gp2, gp3, gp4)
  ggsave("C:/Userdata/Shared/Dofiles/DoAnalysis/JanKalle/ReplicationFiles/AnalysisScripts/Figures/FigureA19.pdf", width = 15, height =10, units = "in") 
  
  modelsummary(Models[c(1, 2, 3, 4)], output = "C:/Userdata/Shared/Dofiles/DoAnalysis/JanKalle/ReplicationFiles/AnalysisScripts/Tables/TableFigureA19.tex")
  
sink()  
  